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Abstract 



rink mode in a periodic 

1 

lj In a line-tied system, 



It is well known that the radial displacement of the m = 1 internal 
screw pinch has a steep jump at the resonant surface where k • B = 0. 
relevant to solar and astrophysical plasmas, the resonant surface is no longer a valid concept. It 
is then of interest to see how line-tying alters the aforementioned result for a periodic system. If 
the line-tied kink also produces a steep gradient, corresponding to a thin current layer, it may 
lead to strong resistive effects even with weak dissipation. Numerical solution of the eigenmode 
equations shows that the fastest growing kink mode in a line-tied system still possesses a jump in 
the radial displacement at the location coincident with the resonant surface of the fastest growing 
mode in the periodic counterpart. However, line- tying thickens the inner layer and slows down the 
growth rate. As the system length L approaches infinity, both the inner layer thickness and the 
growth rate approach the periodic values. In the limit of small e ~ B^jB z , the critical length for 
instability L c ~ e -3 . The relative increase in the inner layer thickness due to line- tying scales as 
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I. INTRODUCTION 



Magnetohydrodynamic (MHD) kink instability in a screw pinch is of great interest and 
has been under intensive study for decades. In a periodic system the internal kink mode is 
unstable if the safety factor q (for a cylinder with axial periodicity length L, q = 27rrB z / LB^) 
drops below unity somewhere within the plasma and increases with radius. Unstable modes 
have helical symmetry, depending on z and </> in the combination kz ± <p, and possess so- 
called resonant surfaces r = r s , on which k-B = kB z ±B^jr s = 0; due to the periodicity, the 
resonant surfaces must also be rational surfaces (if k = 2nn/L, then on r = r s , q = ±l/n 
is a rational number). The driver of the instability lies within r s , and the radial component 
of the plasma displacement nearly vanishes outside it. It has been shown that the thin 
transition layer near r s predicted from linear theory corresponds to an infinite current sheet 
in finite amplitude theory, at least within the framework of reduced, ideal MHD.jlj In a 
plasma with large but finite Lundquist number, the steepening of the current layer must 
trigger resistive energy release, and it has been suggested that this energy release corresponds 
to the sawtooth crash. 

The kink instability has also been proposed clS cL CcLUSG of solar flares. In this scenario, the 
instability occurs in a force free coronal magnetic loop which emerges from the photosphere 
and is twisted by photospheric motions. The possibility of forming a thin current sheet via 
the kink mode is particularly interesting, because the Lundquist number S in coronal loops 
is so high (~ 10 10 ~ 12 ) as to otherwise effectively preclude fast resistive MHD processes. For 
this reason, the instability and its current sheet are also of interest for the coronal heating 
problem. Following Parker's original suggestion, Q] many authors |3,00>E1 have shown that 
random shuffling of the coronal magnetic fieldlines by photospheric motions progressively 
increases the current density, resulting in sporadic energy release. How these spatially 
intermittent currents are produced, and in particular whether MHD instability plays any 
role, remains unclear. The kink instability is an interesting possible model of energy buildup 
and release. Line-tied kinking has also been captured in laboratory experiments; a recent 
example was reported in Ref. 0]. 

However, the theory of the kink mode developed for periodic plasma does not carry 
over to line-tied systems such as coronal loops in a straightforward manner. It is generally 
agreed that line-tying is stabilizing, in the sense that a system with periodic boundary con- 
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ditions is more unstable than an identical system of the same length with line-tied boundary 
conditions. Stabilization arises from the extra tension force associated with the anchored 
footpoints. Roughly, the system should be stable if the Alfven travel time along the cylinder 
is less than the inverse of the growth rate in a periodic system. The first quantitative anal- 
ysis of the effect of line-tying on stability was given by Raadu,|8[ who minimized the energy 
of a restricted class of perturbations. However, there is a deeper issue. Eigenfunctions with 
helical symmetry cannot satisfy boundary conditions on surfaces of constant z. Therefore, 
the significance of rational surfaces in aperiodic systems is open to question. Nevertheless, 
numerical studies have shown that when a line-tied flux tube is twisted 

by motions of its endpoints, it eventually kinks, and the kinked state has a thin current 
layer (whether it is a current singularity or not is very difficult to examine by a numerical 
calculation). On the other hand, it has been shown quite rigorously that current sheets with 
the simple ribbon topology found from kink theory cannot appear in line-tied magnetic fields 
as long as their endpoint motion is a smooth function of position. 14 , 3] Thus, the question 
of what sets the current density in a kinked, line tied system is still open to investigation. 

The purpose of ths paper is to clarify the relationship between the periodic and line-tied 
kink instabilities. In particular, we address the formation of thin current layers in line tied 
systems. In a periodic system, the rational surfaces constitute a natural set of separatrices, 
where current sheets can form. In a line-tied configuration, the notion of rational surfaces 
is absent; as such, it is not clear a priori where a thin current layer or current sheet, if 
any, would form. The present study addresses this important issue. And the answer, in 
short, is that the steep gradient in a line-tied system will appear at the resonant surface 
corresponding to the fastest growing periodic mode in a infinitely long system. A detailed 
comparison between the periodic eigenmode and the line-tied eigenmode is made. Through 
the comparison, the effects of line-tying on the growth rate and the eigenfunction can be 
clarified. 

We limit ourselves to the linearized problem in this study, and leave the nonlinear problem 
to future work. This paper is organized as follows. Sec. |H] give the details of the model 
and the governing linearized equations. Sec. IIHI presents the main results: comparsion of 
the periodic eigenmode and the line-tied one. We mostly concentrate on the fastest growing 
mode in both cases. Details are given for how the latter approach the former in the limit 
when the system length L goes to infinity. An important quantity — the thickness of the 
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internal layer — is found to follow a scaling law in that limiting process. Some other scaling 
laws regarding the critical length, the internal layer thickness and the mode localization 
along the axial direction at marginality are also found. We conclude and give a discussion 
in Sec. II VI Two appendices follow the main text. Appendix [X] gives the details of the 
numerical methods, and a semi-analytic calculation for the critical length and the growth 
rate is detailed in Appendix |Bj 



II. MODEL AND EIGENMODE EQUATIONS 

We assume that the plasma pressure is negligible, p — 0, and that the equilibrium plasma 
density p = const, for simplicity. In cylindrical coordinates (r, <fi, z), the equilibrium magnetic 
field is 

B = B 4 ,(r)4> + B z (r)z, (1) 
which satisfies the force balance equation 

4^_^ = . (2) 
dr \8n J Anr w 

Assuming e 7 * time dependence, the linearized ideal MHD equation can be expressed in terms 
of the Lagrangian displacement £ as 

7 2 P£ = -^(V x (V x (£ x B))) x B + x (V x (£ x B)) = 0. (3) 

It is convenient to decompose the displacement £ into the radial, perpendicular, and parallel 
components as 

£ = £rf + 6,»7 + f[|b, (4) 

where 

f) = (B z 4> - B^)/B, (5) 

and lo is the unit vector along the equilibrium magnetic field. 

In the equilibrium, the <fi direction is ignorable; we assume azimuthal dependence e im ^. 
Taking the b component of Eq. © gives £n = 0. The remaining two independent compo- 
nents of © are: 
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+ij ((v A v A .y - ™±) ?„ ( 6 ) 



(7) 



where primes denote d r ;V^ z — B g f \^np, Va^ — B^j y/Aitp , and Va = B/^/Aup, respectively. 

In this work we are mostly interested in systems with a strong guide field, i.e., Va z 3> Va<$>- 
In that case, the kink mode is nearly incompressible and the growth rate is much smaller 
than the Alfven frequency, i.e., 7 2 <C V^/a 2 , where a is a characteristic perpendicular length 
scale of the flux tube. That means that large terms (~ 0(V^/a 2 ) |£|) on the right-hand 
side (RHS) of Eqs. © and (JJJ) nearly cancel each other; the small unbalanced remainders 
are balanced by the inertia term on the left-hand side (LHS). The cancellation of large 
terms could lead to a significant loss in accuracy when the eigenmode equations are solved 
numerically. To avoid that we reformulate the eigenmode equations as follows. Let 

= + (8) 



mV a V r 

where the first part on the RHS gives (approximate) incompressibility and £ v is a small 
correction. Substituting (JSJ) into the eigenmode equations, we have 

1% = Vlf> + {^f + 2V M V^ £ + (1 "" 2) S r 

\ V mr mrVAz thVaz J J 

irV Az V A4 ,_ d „ + f^%_ + ir^% ~ VLW^ _ ^V^VaA q g 



m ' \mVAz tt^Vaz m 
-V^Vk% - — f + VkV^ J Z &, + 1 6 



+ -(_f^ + M^_^k_M- )£i , ( 9) 
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( ^ ? " ) £»■ ( 10 ) 



In this way the large terms on the RHS have been cancelled explicitly; all terms are 
Oiyl^/a 2 ) |£| or smaller. 

To make the numerical solution easier, we assume the existence of a conducting wall at r = 
r , which is made sufficiently far away that stability is only weakly affected. Two conducting 
plates at z = ±L/2 anchor the magnetic field footpoints. Under these assumptions, the 
boundary conditions are £ r | z=±£ / 2 = £ v \z=±l/2 = 0, and £ r |r=r — 0. The regularity condition 
at r = requires that £ r ~ r H-i anc i £ ^ r \m\+i ag r approaches zero. 



III. NUMERICAL SOLUTION OF EIGENMODE EQUATIONS 

In this section we present the numerical solutions of the eigenmode equations. The 
numerical method is described in Appendix El The model system we solve in this work has 
the following force free equilibrium profiles: 

U^ = erexp(-r 2 /2), (11) 



V Az = v / l + e 2 (l-r 2 )exp(-r 2 ). (12) 

This simple system has only three free parameters: e, L, and tq. We fix the conducting 
outer boundary at r = 5, sufficiently far that it only weakly affects the stability and the 
eigenfunction. The parameter e ~ Va^/Vaz is usually a small parameter in tokamak and 
coronal loop applications; in this work we systematically explore the region e < 1. For a 
given equilibrium, we solve the eigenmode equations for different L. We are particularly 
interested in studying how the growth rate and the thickness of the internal layer vary with 
the system length L. 



A. Periodic Solutions 



Before going into the results for line-tied systems, we present the results of the corre- 
sponding periodic systems. These results will serve as references to the line-tied solutions. 

The analytical theory of the m = 1 ideal internal kink mode in the limit e 1 is given 
by Rosenbluth, Dagazian, and Rutherford in Ref. [l] (hereafter RDR). Key results relevant 
to the present study are summarized as follows: (1) For an unstable eigenmode with spatial 
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dependence exp(ikz + im(fr), with m = ±1, the eigenfunction Cr{ r ) has a steep gradient (the 
internal layer) at the resonant surface r = r s , i.e., where k-B = kB z + mB^/r s = 0. Outside 
of the internal layer, £ r ~ £ a = const for r < r s and £ r ~ for r > r s . (2) The growth rate 
of the unstable mode is approximately (there is a misprint of a factor of two in the original 
paper) 

7T r a 

7 - - J7T, T? \/ 1 „ 3 / 9idr, (13) 



ir. s 







where 

9l = k 2 r(3k 2 r 2 Vl + 2kmrV Az V A<t> - V^). (14) 

(3) The solution within the internal layer is approximately 

1. A 2 _ x fx\(k-V A ) 



e, - 2& ^ - ^ tan 7 J J ' ( 15 ) 

where x = r — r s . 

Some scaling laws can be deduced from these results. First, we have 7 oc e 3 in the 
following sense. In the limit e < 1, V Az is approximately constant. If we vary e and let 
k oc e, the resonant surface will be approximately located at a fixed radius, and the growth 
rate 7 will be proportional to e 3 , since g\ oc e 4 and |(k • "V a)'\ t oc e. Likewise, the thickness of 
the internal layer is proportional to e 2 . To be precise, we define the thickness of the internal 
layer A as the distance between the two radii where £ r = (3/4)£ and where £ r = (l/4)£ . 
Eq. (fT5j) gives 

A^2 7 /|(k-V A )'| rs . (16) 

Hereafter we will use the subscript "0" to denote properties of or related to the fastest 
growing periodic mode: for a given equilibrium, the fastest growing mode appears at the 
wavenumber k = k$, with the growth rate 7 = 70, and the internal layer thickness Ao- From 
the scaling laws we have k oc e , 70 oc e 3 , and A oc e 2 . 

Figure shows the growth rate 7 as a function of k, for different e. Solid lines denote 
the growth rate calculated by the code, and dashed lines are the approximate growth rate 
calculated from Eq. ([13)1 . As expected, the agreement becomes better for smaller e. For 
each e, the growth rate peaks at k — ko, and the corresponding 70 follows the scaling law, 
70 oc \ko\ , as indicated by the dashdot line. The numerical values of ko, 70 , and A for 
different e are summarized in Table HI Figure El shows the eigenfunctions of the fastest 
growing modes for e = 1, 0.5, 0.25. The radial displacement £ r of each shows a jump at 
r ~ 1. The jump becomes steeper, and the twist (C, v ) becomes more localized, for smaller e. 



Notice that here the axial wavenumber k is treated as a continuous parameter. This 
treatment actually corresponds to periodic eigenmodes in an infinitely long system. For a 
periodic system with a finite length L, the wavenumber k is quantized as k = 2mr/L, with 
integer n. As such, only a finite number of unstable modes are present in a finite length 
system. 

B. Line-tied Solutions 

We start by observing some general characteristics of the eigenmodes in line-tied systems, 
then move on to the scaling laws obtained from analyzing the numerical solutions. We focus 
only on the fastest growing mode for a given configuration in the first two sections. Higher 
harmonics are briefly discussed in section IIIIB3I 

1. General Observations 

Figure El shows the eigenfunctions of the fastest growing mode for e = 0.5, L = 300, 
which demonstrate some general characteristics observed throughout all the cases we have 
tried. The first thing to notice is that there are many oscillations along the z direction. The 
wavelength is approximately the same as the corresponding fastest growing mode in the 
periodic case. Second, the radial displacement also has a jump at r ~ 1, the same location 
as the jump of the fastest growing periodic mode. Third, the eigenmode is more or less 
localized to the center, rather than being broadened out to the whole z domain, with an 
envelope ~ cos(7r;z/L), which seems to minimize the field line bending due to line-tying. In 
fact, the cos(tcz/L) envelope is what Raadu used in his energy principle analysis, [8] yet the 
numerical solutions suggest otherwise. We will come back to this issue later. 

Since the wavelength along z is approximately the same as that of the fastest growing 
periodic mode, one may "filter out" the fast oscillations along z by dividing the solutions 
by exp(ikoz). The results are shown in Figure 0] The remaining "envelopes" of the eigen- 
functions become slowly varying along the z direction. This feature has been utilized into 
the choice of the basis functions of the numerical method, detailed in Appendix El in short, 
instead of having to resolve the fast oscillations along z, we only need to resolve the slowly 
varying envelopes, therefore many fewer basis functions are needed. 
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Figure 03 shows the fastest growing mode for different L, with a fixed e = 0.5. We 
observe that as L becomes larger, the jump of the radial displacement becomes steeper, and 
the eigenfunction becomes broader along z. For a really long L, the envelope cos(nz/L) 
becomes a good approximation. Figure El shows the eigenf unctions at z — 0, for various 
L, as compared to the periodic fastest growing mode. As L becomes larger, the line-tied 
eigenfunctions at the midplane approach the periodic ones. For the case L = 1500, which is 
not shown on the plot, the eigenfunctions at the midplane are virtually indistinguishable from 
the periodic ones. The parallel component of the perturbed current Jii = b - 
at the midplane is shown in Figure We see the thin current layer of the periodic solution 
is smoothed in line-tied cases. As the system length becomes longer, the line-tied solution 
approaches the periodic one. It should be pointed out, however, that the aspect ratio of 
these systems are much larger than in natural systems such as coronal loops. 

Additional insight may be obtained by decomposing the eigenfunctions into Fourier har- 
monics. Because the wavelength along z is approximately the same as the wavelength 
corresponds to k , we consider the following "shifted" Fourier decomposition: 

oo 
n=— oo 

Figure|Hl(a) shows the six most significant Fourier components for the case e = 0.5, L = 300. 
We observe that they roughly form three pairs. The two components in each pair have 
approximately the same amplitude such that they nearly cancel each other at the ends, 
z = ±L/2. However, the expanded view (b) about the jump at r ~ 1 shows that each 
component has a jump at its own k • B = surface, precluding cancellation within the 
narrow internal layer. As a result, many (actually an infinite number of) Fourier harmonics 
are needed to achieve a full cancellation at the ends. For a longer system, the distance 
between neighboring k • B = surfaces becomes smaller. As shown in Figure |H1 (c)(d), for 
the case e = 0.5, L = 1500, the two dominant components are nearly identical throughout the 
whole domain, and the components of other harmonics are very small. One may anticipate 
that in the limit L — > oo the solution is completely dominated by two components, which 
give a z dependence ~ exp(zfenz) cosjirz / L). This is similar to the external kink mode 
solution calculated by Hegna lg and Ryutov et. al. 
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2. Scaling Laws and "Universalities" 

It would be desirable if we could deduce some scaling laws from the numerical solutions. 
In particular, we are interested in how the internal layer thickness, the growth rate, the 
mode localization along z, the critical length for instability, etc., vary with e or L. To be 
precise, let us define the quantities to be measured. The thickness A of the internal layer is 
measured at z = 0. Let £ a = |£ r | r=0z=0 . The thickness A is defined as the distance between 
the radii where |£ r | = (3/4)£ and |£ r | = (l/4)£ . The length I of an eigenmode along z is 
measured at r = 0, and defined as the distance between the two points where |£ r | = (l/2)£ a . 
We use the subscript "c" to denote properties associated with marginal stability: L c is the 
critical length for a given e; l c is the length, and A c is the internal layer thickness, of the 
marginally stable eigenmode. 

Table HP summarizes L c , A c , and l c for various e. For small e, we observe the following 
scaling laws: L c oc 1/e 3 , (A c — A )/A oc 1/e, and l c j L c oc e. The scaling L c oc 1/e 3 can 
be understood as follows. From a semi-analytic calculation (see Appendix El for details), we 
have an estimate for the critical length 

L c ~2\V Ag \ r J>y (18) 

where Va z is evaluated at the resonant surface of the fastest growing periodic mode. For 
e =1, 0.5, 0.3, 0.25, 0.2, 0.15 this estimate gives L c -22.9, 130, 538, 910, 1746, 4079. 
Compared with the L c summarized in Table |Hj the agreement is very good for small e. 
Since Va z — const and 70 oc e 3 for small e, the scaling law follows. In the astrophysics 
literature, more often the critical twist of the magnetic field for a given flux tube length is 
considered. Let us consider the twist at the critical length: 

A0 C oc ^± oc i oc Ll'\ (19) 
rB z e A 

If we consider the flux tube length L as given, and twist up the magnetic field to make it 
unstable, the critical twist scales as L 2 / 3 . 

The quantity (A c — A )/A measures the relative increase in the internal layer thickness 
due to line-tying, and l c /L c measures the localization along z of the marginal mode. The 
scaling law (A c — Ao)/Ao oc 1/e indicates that the increase in internal layer thickness is 
greater for smaller e, at marginal stability. And the scaling law l c /L c oc e indicates that 
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the marginally stable mode becomes more and more localized along z for smaller e. Figure 
shows |£ r | of the marginally stable eigenmodes for various e. The localization along z is 
evident for smaller e. Note that the trend of localization along z and broadening of the 
internal layer are reciprocal to each other. This is consistent with the understanding based 
on Fourier mode decomposition — localization along z means a lot of Fourier harmonics are 
involved, that also means significant broadening of the internal layer. 

Let us now look at the internal layer thickness A as a function of L. We have already 
observed that A — > A as L — > oo. In an attempt to reveal some "universality", we plot the 
functional dependence for various e, with a proper choice of variables and normalizations, on 
the same diagram. A natural choice of normalization for L is with respect to L c ; therefore 
we choose L c / L as the horizontal axis. The choice for the vertical axis is not all that obvious. 
Since (A c — A )/A oc 1/e, we choose e(A c — A )/A as the vertical axis. With this choice 
of variables, the marginally stable solutions for different e all appear roughly at the same 
point. Figure ITUl shows the resulting plot in log — log scale. We observe that all curves for 
different e roughly coincide, and follow the scaling law e(A — A )/A oc (L c /L) 2 - 5 . The 
curves deviate from the scaling law near marginality. 

Finally, let us look at how the growth rate varies with L. It is found that 7 — > 70 as 
L — > 00. We again try to plot the relationship between 7 and L through a proper choice 
of variables. A natural choice is L c /L as the horizontal axis, and 7/70 as the vertical axis. 
This choice brings all marginally stable modes to the same point (1,0) on the plot. Figure 
mi shows the resulting plot, for various e. Once again, all curves roughly coincide. The 
relation between 7/70 and L c /L can be obtained by the semi-analytic calculation detailed 
in Appendix iBl as 

L c _ 7/70 (2Q) 



L tanh 1 (7/7o)' 

which is in a good agreement with the numerical result. The line-tied growth rate 7 ap- 
proaches the fastest periodic growth rate 70 rather quickly after the critical length is ex- 
ceeded. For L = 2L C , about 95% of the periodic growth rate is recovered. 

3. Higher Harmonics 

Thus far we have been focused on the fastest growing mode. Here we briefly discuss the 
higher harmonics. Figure El shows four harmonics of the case e = 0.5, L = 200. Since 
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the length L is not much longer than the critical length L c = 142.4, only three unstable 
modes are present. The kinking central column breaks into blobs for higher harmonics; 
the growth rate decreases with the increase in the number of blobs (hence more field line 
bending). Figure El shows the four most unstable modes for the system with the same 
e = 0.5, but much longer L = 1000. The higher harmonics still show the same general 
characteristics, but the growth rates of them are very close to the fastest one. In this case 
the physical significance of the fastest growing eigenmode becomes questionable, since those 
higher harmonics are likely to play an equally important role when the instability develops 
and reaches nonlinear saturation. 



IV. SUMMARY AND DISCUSSION 



In this paper we address the issue of the internal kink mode in a line-tied system, as 
compared to that in a periodic system. We find that the fastest growing internal kink mode 
in a line-tied system possesses a steep internal layer, as does its periodic counterpart. The 
internal layer locates at the resonant surface of the fastest growing periodic mode in an 
infinitely long system. Therefore, even though the notion of rational surfaces is absent in 
a line-tied system, we still have a rule of thumb as to where the current sheet would be. 
Line-tying decreases the growth rate and increases the internal layer thickness; only in the 
limit L goes to infinity does the line-tied mode approach the periodic one. 

In the small e limit, the critical length can be estimated rather accurately using Eq. 
flUJ), an d the dependence of the growth rate on L is accurately represented by Eq. (j2*Uj) . 
Dimensionally, the critical length agrees with the physical picture that the system becomes 
unstable when the Alfven travel time along the cylinder is longer than the inverse of the 
growth rate in a periodic system. A similar estimate can also be derived by applying the 
RDR asymptotic analysis to Raadu's energy principle calculation, which gives 

L c * (21) 

7o 

Functionally this is similar to Eq. ()18j) . but significantly overestimates the critical length 
by a factor of tt/2. The reason is that Raadu assumed a cos(n z/L) envelope in his trial 
function, which is not a good approximation for the marginally stable mode, as the latter is 
highly localized to the center. 
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Several interesting scaling laws are deduced from the numerical solutions, most notably 
that L c oc 1/e 3 , l c oc 1/e 2 , (A c - A )/A oc 1/e, and e(A - A )/A oc (L c /L) 2 - 5 . However, 
the computational overhead increases rather rapidly as one decreases e, due to the stringent 
scalings 70 oc e 3 and A oc e 2 . As a result, the attainable parameter range of the numerical 
solution is limited. It would be desirable if the problem could be solved asymptotically in 
the limit e < 1, thereby the scaling laws can be "explained". Our attempt at asymptotic 
solution only succeeds partially. We hope that the numerically found scaling laws can give 
some hint as to how the full analysis could be carried out. 

As for application to realistic situations, Eqs. fT£j) . (J2U|) . and the observed scaling laws 
are only valid in the small e limit. For small e, the critical aspect ratio L/a can easily go 
up to hundreds or thousands in order to make the flux tube unstable, and one may have 
trouble finding such a long flux tube in nature. Furthermore, for a naturally arising flux 
tube twisted up by footpoint motions, one can anticipate that the tube length will not be 
much longer than the critical length (or equivalently the twist will not be much greater then 
the critical one) before the instability sets in and completely alters the system. Therefore, 
even though the line-tied eigenfunction does approach the periodic one as L — > 00, the 
physical significance is not clear. All these considerations point to the much more important 
nonlinear problem. Our next step will be to look at the nonlinear equilibrium after the 
instability, which will serve as a direct comparison to the nonlinear equilibrium given by 
RDR. Presently, there is still no consensus as to whether an infinitely thin current sheet or 
a current layer with finite thickness would form in a line-tied configuration. And we hope 
we can further clarify this issue. 
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APPENDIX A: NUMERICAL METHOD 



We discretize the eigenmode equations by a Petrov-Galerkin scheme |18j|. For a generalized 
eigenvalue problem of the form XCiu = C2U, the weak form of the same system is \(C\U, w) = 
(C2U, w) for all w in a suitable space of test functions, where (•, •) denotes the inner product 
over the spatial domain. For our numerical method, we will approximate u and w with finite 
dimensional spaces, therefore turning the original partial differential eigenvalue equations 
into a matrix eigenvalue problem. In particular, we expand the variables £ r and as 



2N r N z 

EX' 



ip? (r) 



(Al) 



Substituting this expression into the eigenmode equations, we then require the resulting 
equations to be satisfied identically when projected onto the 2N r x N z basis functions of the 
test space 



<p*(z). 



(A2) 



For our convenience, we take 







r 










(«- + 9 ' 







P(z) 



(A3) 



(A4) 



and the inner product is defined as 



h 

9i 



h 

92 



pro 

dz / dr(f 1 f; + g^*). 
Jo 



L/2 r ro 

-L/2 



(A5) 



The choice of ipr and gives the correct form of energy integration in the r direction, 
and the factor 1, 



'1 — (2z/L) 2 in (j) u allows us to utilize the orthogonality of Chebyshev 
polynomials, which are used to construct the basis functions (see details below). 

After lengthy calculations and integration by parts, the discretized approximation to the 
eigenmode equation can be written as 



Ufa, 



(A6) 
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-rVl + 



A 7 = 



r 2 V A V Az 



m 



A, 



rV A V Az 



m 



A 9 = rV\- 
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TQ 



dr dr 



dr n dr 



dr. 



(A18) 



, 2r 2 V Az V M 
A io — • 



m 



(A19) 



A 



11 



2(m 2 - l)V M V Az 
m 



(A20) 



A 12 = -rV A(j) V A , 



(A21) 



A i3 = V A V A4> ; 



(A22) 



pro 

Mf = 

Jo 



d^ dip? 



(A23) 



(A24) 



A 15 = 



(A25) 



A m ^ 

^-16 — • 



(A26) 



L/2 J2iv 
■L/2 dz 



(A27) 



f L/2 T'Wj 

J — _, 



L/2 
L/2 

t 

L/2 



dz 



^ <b v dz. 



(A28) 



(A29) 
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We have tried various basis functions. A straightforward choice satisfying boundary 
conditions and regularity conditions at r = based on Chebyshev polynomials would be 



( r ) 

V>£(r) 



( r / ro )H-i[ T/i ( 2r / ro _i)_i] 








( r / ro )|m|+i T _ iVr _ 1 (2r/r - 1) 



when fj, < N r 



when fi > N r , 



and 



P{z)=T u+1 (2z/L)-T„- 1 (2z/L), 



(A30) 



(A31) 



(A32) 



where T n denotes the Chebyshev polynomial of order n.[19] This choice of basis functions 
works well, but a significant amount of basis functions are needed, as the eigenmode typically 
possesses a steep gradient at some radius r = s and many oscillations along z. To remedy the 
steep gradient problem, we need to pack more resolution about r = s. This is implemented 
by successively applying three arctan/tan mappings, which are recommended in Ref. jl^ . 
as follows: First we map r G [0, ro] — » r\ G [0, 1] with s — ► 1/2 by 



7T 



tan ai tan 



7r r 
2^1 



with 



tan 



2 r 



(A33) 



(A34) 



Then a second mapping from [0, 1] — > [0, 1] 

1 



r 2 



tan 



a2 tan hr ri - 



1 

2 



(A35) 



7T v V V 2 y 

is applied, which packs more resolution about r x = 1/2 (which is mapped into r 2 = 1/2), 
with larger a 2 for more packing. Finally, we map r 2 G [0, 1] — > r% G [0, 1] with 1/2 — > s' by 



r 3 



2 

— tan 

7T 



-1 



a 3 tan ( -r 



7T 



(A36) 



with 



a 3 = tan . (A37) 

In this mapping, s' can be any number between and 1, but we usually choose some s' > 1/2. 
This is motivated by the observation that typically the eigenmode structure is rather trivial 
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(~ 0) in the region r > s, therefore one can put less resolution there. The basis functions 
are then defined in terms of r 3 as 



and 







y*. 




Vr 




y%_ 





[r 3 



iH-i 



[T M (2r 3 - 1) - 1] 




when fi < N r 



(A38) 







when fi > N r . 



(A39) 



(r 3 )H+ir M _ JVr ._ 1 (2r 3 -l) 
To resolve the many oscillations along z without too many basis functions, we use the 
observation that the fastest growing mode, upon dividing by e tk ° z with ko being the wave 
number of the fastest growing mode in the corresponding infinite system, yields a smooth 
function with length scales ~ L. This suggests that one may use 



<j> v (z) = e lk ° z (T u+1 (2z/L) - T v ^{2z/L)) 



(A40) 



as the basis functions. This proves to be very efficient. 

The integration along r is carried out numerically by Gaussian quadrature. The integra- 
tion along z can be done analytically, which gives 



D1 



7T 



3/2 
1 
-1/2 

-1/2 



1 

-2 2 
-3 3 
-4 



1/2 

-1/2 

1 



(A41) 



(A42) 
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or = ^ x 



-2n(ji + l) 
-4/i 



\i = v 
H < v>, v — ji even 



+ 2ik n 



1 
-2 2 
-3 3 
-4 



1,2 n /J.is 



(A43) 



The resulting finite dimensional eigenvalue problem is then solved by a Jacobi-Davidson 



21. 



iteration method, which is described in detail in Ref. 

The growth rates of unstable modes are benchmarked with an initial-value MHD code, 
M MHOD;0 and the critica! lengths for in3tabiH ty are a ls o in g ood a„t with the 
results published in Ref. Q|. A periodic version of the code is also developed, in which 
spatio-temporal dependence of the form exp(ikz + imcf) + -yt) is assumed. 

APPENDIX B: SEMI- ANALYTIC ANALYSIS OF LINE-TIED INTERNAL 
KINK 



The linear analysis of the periodic internal kink modes in RDR goes as follows. First the 
eigenvalue problem is approximated by a single equation for £ r . The equation is then solved 
approximately in the internal layer and the outer region separately. Finally, asymptotic 
matching of the inner and the outer solutions gives the growth rate of the unstable mode. 
In this appendix, we try to extend the analysis to a line-tied system. An estimate for the 
critical length and the relation between the growth rate and the system length are derived. 
It should be pointed out though, that the analysis is not completely satisfactory. At certain 
point, we are forced to appeal to numerical work to validate certain steps. That is why we 
call it semi- analytic. 

Following the steps of the analysis of RDR, first we have to derive a single equation for £ r . 



We start with the variational principle for the eigenvalue problem. 



251 ] Assuming m = ±1, 
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for a displacement £, the energy variation can be written in term of £ r and £ v as 

/ L / 2 r r ° / * \ 

dz / dr r 3 \V A ■ V£| 2 + A 6 \d z £ r \ 2 - iA 7 (d^ d& - c.c.) 
- L /2 J 



iA 8 ( d z £„ d z £ r - c.c. ) + A. 



-A 12 (L*d z C + c.c.) - Ats [L*d z ^ r + c.c.) - A 







2 


)-A 16 


in 





:bi) 



and the inertia term can be written as 

•L/2 



dz drr A 1 \C\ 2 + A 2 |£ r | 2 - iA 3 (£, C ~ c.c. 

L/2 Jo 



-iA A [L £ r - c.c. + 







2 


)+ 


in 





(B2) 



where the coefficients A\ — A 16 are defined in Appendix |XJ and c.c. denotes complex conju 
gate. The eigenvalue problem is equivalent to the following variational problem: [2 

,2 W(t) 



7 



and 



(B3) 



(B4) 



Since £ v is small for kink modes in the small e limit, we may approximate as 

-L/2 r r 



3 \C'\ 2 



dz / drr 

-L/2 JO 



(B5) 



That is, that only weakly depends on ^ ; therefore we may minimize W($,) with respect 
to £ n first in the variational problem ()B4|) . The minimization gives the following relation 
between £ r and £ v : 



(A 9 d 2 z + A 16 )Z V = iA 7 d% + iA & d% - A 12 d z C - A l3 d z £ r . 



(B6) 



Since rd z ~ ikr ~ e, this shows that £ v ~ e 2 ; therefore the approximation ()B5j) is self- 
consistent. Using Eq. (|B6|) in (|B3|) to eliminate as many £ v as possible, we may write W 

as 

-L/2 



ro 



efc / dr [r 3 \V A - Vi'X + As\d z i T 

L/2 JO 



A, 



+ A 



16 



6= 



• (B7) 
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To complete the minimization, we still have to solve ()B6j) for £ v in terms of £ r . This could 
formally be done in terms of the Green's function, but let us proceed with the following 



approximation. We expand W up to 0(e 4 ). The term Ac 



in the integrand is 0(e 6 ) 



therefore is negligible; and from Eq. ()B6 



i v = (iA 7 dlC + iA 8 d% - A 12 d z C - A 13 d z £ r ) /A 16 + 0(e 4 



(B8) 



After some algebra and integration by parts in r, we have 

•L/2 



W 



pro 

dz I dr 

L/2 JO 



r 3 \Va ■ Vei 2 - r 5 \d z {Y A • VQ\ 2 + 3rVj 2 I^J 2 - rV" \d z £ r \ 2 



+imV Az V A<t> r 2 (rd z Cd% - rd z C r d% + ^CUr ~ c.c) 



+ 0(e 6 ). 



(B9) 



Integration by parts in z gives 

"L/2 pro 



dz / dr imV Az V A(j) r 2 {rd z £'*d% - rd z Cd 2 z C - c.c.) 

-L/2 JO 

= / drtmV Az V M r 2 (rd z Cd z Z r -c.c.) L/2 . (BIO) 



This term may be negligible in Eq. (jB9|) . as one can estimate 

(■L/2 

(rd z C*d z £ r - c.c.fJL : / dz (d 2 z £d z £ r - c.c.) ~ 1 : kL ~ e 2 : 1, (Bll) 

J -L/2 

where in the last step we use the observed scaling law that to have unstable mode, L > 
L c oc 1/e 3 . We may also neglect the second term in the integrand of Eq. (jB9|) . since 



r 5 \d z (V A ■ VC)| 2 : r 3 \V A ■ V£| z ~ k 2 r 2 : 1 ~ e 2 : 1. 



/|2 . ,2 2 



(B12) 



Therefore W can be further simplified to 

«L/2 



w 



pro 
dz / dr 

L/2 Jo 



r 3 \V A ■ V£f + 3r 3 V 2 z \d%\ 2 - rV 2 A(f> \d z £ r 



+imV Az V M r 2 {d 2 £dz£r - c.c.) 



(B13) 



Now we have the expressions for the energy W and the inertia / in terms of £ r alone, Eqs. 



()B13|1 and ()B5j) . We may now substitute the expressions into the variational principle (|B4|) to 
obtain an single partial differential equation which approximates the full eigenvalue problem. 
There is a subtle point here. The form of the approximate W requires four boundary 
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conditions for £ r in the z direction. Two apparent ones are the line-tied boundary conditions 
£ r \ z =±L/2 = 0. A possible choice of another two boundary conditions is d z ^ r \ z= ±L/2 = 0, 
which makes all boundary terms vanish when turning Eq. (|B4J) into a partial differential 
equation through integration by parts. Assuming that, the approximate eigenvalue equation 
is 

(r 3 (- 7 2 + (V A ■ V) 2 )£)' + (3rVjA 4 + 2imr 2 V Az V A ^ z + rV 2 M d 2 z ) £ r = 0. (B14) 

It is hard to justify this particular choice of boundary conditions rigorously. In particular, it 
is not the same as requiring £ v = at the boundary (see Eq. (|B8J) ). Nonetheless, numerical 
solutions show that this new eigenvalue problem is indeed a good approximation to the 
original one in the small e limit. As we will see, it turns out that £ r is vanishingly small 
within the boundary layers near z = ±L/2, therefore indeed d z ^ r \ z= ±L/2 — 0. 

One nice feature of Eq. (jB14|) is that, if we consider a periodic problem by letting d z — ► ik, 
the approximate equation of RDR is recovered. This allows us to apply a similar asymptotic 
technique to the present problem. Motivated by the numerical solution, let us write the 
eigenfunction as £ r (r,z) = h(r, z)e lkl)Z , where k is the wavenumber of the fastest growing 
periodic mode, e tk ° z is the fast oscillating part, and h(r, z) is the envelope slowly varying 
along the z direction. Substituting this into the energy Eq. (|B13j) . the dominant term in the 
integrand is the positive definite r 3 \V a • V^| 2 = r 3 \ {VA z d z + zk • Va) h'\ 2 ~ 0(e 2 ), and all 
other terms are 0(e 4 ). For an unstable mode, which requires W < 0, we must have h' ~ 
everywhere except where ko • ~ 0. Therefore, the internal layer occurs at r = r s , where 
ko • = 0; outside of the internal layer, h(r, z) ~ h (z) when r < r s , and h(r, z) ~ when 
r > r s . 

Integrating Eq. ()B14|) once along r, we have 

(-7 2 + (Va ■ V) 2 ) e r (r,z) = 

_ J_ f df (3f 3 Vl(f)dt + 2imf 2 V Az {f)VA<t>{f)dl + fV^d 2 ) ^(f, z). (B15) 
r Jo 

Using £ r (r, z) = h(r, z)e tk ° z in Eq. (|B15|I . in the region r < r s we have 

(-7 2 + (VaA + 2k • V A ) 2 ) h! ~ ^4^- f gi (f)df, (B16) 

r Jo 

where 

9l (r) = 3k 4 r 3 V 2 z + 2mkyVA z VA 4> - fcgrV^. (B17) 
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In the RHS of Eq. (|B16|) we make use of that h(r,z) ~ ho(z) when r < r s , and all 
derivatives on the envelope function are neglected as they are smaller compared to deriva- 
tives on the fast oscillating part e lk ° z . Within the internal layer, we may assume that 
(— 7 2 + (VA z d z + ik ■ V^) 2 ) h! is nearly independent of r; that is 

(-7 2 + (V Az d z + tk • V A ) 2 ) b! = f(z). (B18) 

The function f(z) can be determined from the outer region by evaluating Eq. (|B16|) at 

r = r s : 

f(z) = ^P- n 9l (r)dr. (B19) 

r s JO 

Furthermore, Vaz can be approximated by Vaz(t s ), and k • ~ (k • Y a)'\ Ts x ~ 
— VAz( r s)('{r s )x, where x = r — r s and 

«') ^ < B2 °> 

Using these approximations, the governing equation within the internal layer is 

{yL{d a -igx)*->f)h' = f{z). (B21) 

Here Va z and C' are evaluated at r = r s . Eq. (|B21|) can be formally solved using the Green's 
function: 

_ e if *(*-*) sinh il{z> _ L/2)/V Az ) sinh ( 7 ( g< + L/2)/V Az ) 
l ' ' ] ~ lV Az sinh( 7 L/^) ' ^ Zl) 

where z > {z<) is the larger (smaller) one of z and z. The solution is 

fL/2 

h'(x,z)= G(x;z,z)f(z)dz (B23) 

J-L/2 

The solution in the inner region has to match the solution in the outer region asymptotically. 
That requires that the jump in the inner solution as x goes from — oo to oo be equal to the 
jump in the outer solution from the region r < r s to the region r > r s : 

dxh'(x,z) = -h Q (z). (B24) 

The only x dependence in h'(x, z) comes from the factor ^ x i z -^) in the Green's function, 
and the integration can be done easily: 

°° - 2-7T 

e Kx(z-z) dx = « S ( z _ ^_ ( B25 ) 



\c 
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Using Eqs. JB25J), (lB23|) . (p22|) . and (pT9j) in (1B24L the asymptotic matching condition is 

2?r sinh ( 7 (z - L/2)/V^) sinh ( 7 (z + L/2)/V^) £ gi(r)dr _ 

T77T t7 ■ u r /t/ a 3 ~~ n o{ z ■ [iiZb 

|C | 7Va 2 sinh 

It seems we are now in trouble. If we eliminate ho(z) from both sides, we end up requiring 

that a function of z equals to unity, which is impossible. However, notice that in the limit 

jL/Vaz 3> 1, we have 

sinh ( 7 (z - L/2)/V Az ) sinh ( 7 (z + L/2)/V Az ) 1 



7 Vk z sinh {iL/V M ) 2 7 |Va*| (B2?) 

for most of the domain z G [—L/2, L/2], except in the two boundary layers of thickness ~ 

Vaz/i near both ends, where it drops to zero quickly. Therefore, as long as h (z) is confined 

to the central region outside of the boundary layers, asymptotic matching is possible. Indeed, 

this qualitatively agrees with the observation on numerical solutions, exemplified in Fig. El 

for the case e = 0.5 with different L. On the other hand, at marginal stability, 7 — > 0, and 

sinh ( 7 (z - L/2)/V Az ) sinh ( 7 (z + L/2)/V Az ) {z - L/2) {z + L/2) 

jV Az sinh(jL/V Az ) ~* LVl ' 1 ] 

In this limit, the "boundary layers" cover the whole domain, and we conclude that the 

marginal mode must be highly localized to the center. This is also in agreement with the 

numerical solutions. 

Based on these observations, it may not be unreasonable to eliminate Hq{z) from both 
side of Eq. f|B26|l . then approximate the left-hand side by its value at z — 0. After some 
algebra, that gives the relation between 7 and L: 

vr tanh( 7 L/2|V^| r J / rs g x (r)dr 1 



\C'VazI 7 



3 



(B29) 



Let us consider some limiting cases of Eq. ()B29j) . In the limit L — > 00 



7T 



7 = -\CV^J (B30) 

This essentially recovers the approximate periodic growth rate (|13|) since (k • V^)'^ ~ 
—V Az (r s )('(r s ). For marginal stability, let 7 — > 0, we get an expression for the critical length 

r 2 I^U 



7o 

Finally, we can rewrite Eq. ()B29|l in terms of 7 o and L c : 



(B31) 



^tanhf^ ) = 1. (B32) 
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In this appendix, we get good approximations for the critical length L c and the relation 
between 7 and L without knowing the eigenmode structure along z. We also have no 
information about the thickness of the internal layer. This is not too surprising, as it is 
well known that getting a good approximation for the eigenfunction is more difficult than 
getting a good approximation for the eigenvalue. This can be understood from the variational 
principle, Eq. (jB4|) : Any deviation of the trial function from the true eigenfunction only 
leads to an error of second order in the eigenvalue. We learn from numerical solutions that 
the mode structure along z and the internal layer thickness are correlated — the more a 
mode is localized in z, the thicker the internal layer is — therefore they must be solved all 
together. Attempts on extending the present analysis to resolve the mode structure have 
been unsuccessful. We hope this appendix can serve as a starting point for a more thorough 
analysis in the future. 
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\ k o\ A 


To A 3 


AoA 2 
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1 
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0.124 
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0.03152 


0.2 


-0.1337 
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1.177 x 10" 3 


0.02005 
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-0.1009 


4.908 x 10" 4 


0.01111 


0.6727 
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0.4938 


4.978 x 10" 4 


0.01122 
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e 


L c 


A c 


lc 


lc/L c 


(A c - A )/A 


L c e 3 


e(A c - A )/A 


l c /L c e 


1 


31.17 


0.5605 


15.64 


0.502 


0.4342 


31.17 


0.434 


0.502 


0.5 


142.4 


0.3304 


44.97 


0.316 


1.715 


17.8 


0.858 


0.632 


0.3 


556.8 


0.1842 


115.9 


0.208 


3.132 


15.03 


0.940 


0.694 


0.25 


932.2 


0.1502 


164.8 


0.177 


3.859 


14.57 


0.965 


0.707 


0.2 


1773 


0.1153 


262.4 


0.148 


4.832 


14.18 


0.966 


0.740 


0.15 


4116 


0.08371 


468.5 


0.114 


6.535 


13.89 


0.980 


0.759 



TABLE II: Summary of L c , A c , and l c for various e, as well as the scaling laws they follow. 




-k 



FIG. 1: The growth rate 7 as a function of k, for different e. Solid lines are the growth rates 
calculated by the code, and dashed lines are the approximate growth rates calculated from Eq. 
(|13|). As expected, the agreement becomes better for smaller e. The dashdot line indicates the 
scaling law, 70 oc |/cq| . 
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FIG. 2: Periodic eigenfunctions of the fastest growing modes for e = 1, 0.5, 0.25. The radial 
displacement £ r has a jump at r ~ 1. The jump becomes steeper, and the twist (^) becomes more 
localized, for smaller e. Note that only the solutions within < r < 2.5 are shown; outside of that 
the solutions are vanishingly small. The eigenfunctions are normalized such that £ r (0) = 1. 
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FIG. 3: The fastest growing mode for the case e = 0.5, L = 300. The two panels on the left show 
the real part of £ r and i£ v . The imaginary parts are not shown; they are similar to the real parts, 
only with different phases. The two panels on the right show the module of each component. The 
eigenfunctions are normalized such that £ r | r=z=0 = 1. 
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Re(^ r /exp(ik z)) 



lm(^ r /exp(ik z)) 




FIG. 4: The fastest growing mode for e = 0.5, L = 300, divided by exp(ikoz) to filter out the fast 
oscillations. The remaining "envelopes" are slowly varying functions along the z direction. 
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ig, E=0.5, L=150 



e=0.5, L=300 




FIG. 5: Contour plots of |£ r |for e = 0.5, L = 150, 300, 700, and 1500. All solutions are normalized 
to £ r | r=2=0 = 1. The contours correspond to |£ r | = 0.1, 0.2, 0.8. As the system length becomes 
longer, the jump in |£ r | becomes narrower, and the eigenfunction becomes broader along the z 
direction. For a really long system, the cos(7rz/L) envelope becomes a good approximation. 
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FIG. 6: Eigenfunctions at the midplane, z = 0, for different L, as compared to the fastest growing 
periodic eigenfunctions. e = 0.5. 
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FIG. 7: The perturbed parallel current at the midplane z = 0. The thin current layer of the 
periodic case is smoothed in line-tied cases. The line-tied case approaches the periodic one as L 
becomes larger. 
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FIG. 8: (a) The six most significant "shifted" Fourier components for the case e = 0.5, L = 300. 
(b) Expanded view of the jump at r ~ 1. Vertical lines indicate the k-B = surfaces corresponding 
to each k. (c) and (d) are the same as (a) and (b), for a longer system with L = 1500. 
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FIG. 9: \£ r \ of the marginally stable eigenmodes for e = 0.5, 0.3, 0.25, 0.2. The eigenfunctions are 
normalized to £ r | r =o,z=o = 1- The contours correspond to |£ r | =0.1, 0.2, 0.8. 
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FIG. 10: e(A — Aq)/Ao as a function of L c /L, for various e, in log — log scale. 
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FIG. 11: 7/70 as a function of L c /L, for various e. The dashed curve is analytic relation, Eq. (|2U[). 




FIG. 12: |£ r | of four harmonics for the case e = 0.5, L = 200. Only three of them are unstable. 
The critical length is L = 142.4. 
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FIG. 13: |£ r | of four unstable modes for the case e = 0.5, L = 1000. 
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FIG. 14: Envelope of £ r at r = for the case e = 0.5 with various L. The envelope function is 
confined to the central region and is vanishingly small within the boundary layers near the ends. 
The envelope function approaches cos(tt2;/L), which is shown as a reference, in the limit L — > oo. 
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